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Resonant tunneling through quantum dot under a finite bias voltage at zero tempera- 
ture is investigated by using the adaptive time-dependent density matrix renormalization 
group(TdDMRG) method. Quantum dot is modeled by the Anderson Hamiltonian with the 
1-D nearest-neighbor tight-binding leads. Initially the ground state wave function is calculated 
with the usual DMRG method. Then the time evolution of the wave function due to the slowly 
changing bieis voltage between the two leads is calculated by using the TdDMRG technique. 
Even though the system size is finite, the expectation values of current operator show steady- 
like behavior for a finite time interval, in which the system is expected to resemble the real 
nonequilibrium steady state of the infinitely long system. We show that from the time intervals 
one can obtain quantitatively correct results for differential conductance in a wide range of 
bias voltage. Finally we observe an anomalous behavior in the expectation value of the double 
occupation operator at the dot {n^ui) as a function of bias voltage. 

KEYWORDS: quantum dot, nonequilibrium, Kendo effect, adaptive time-dependent DMRG 



1. Introduction 

Recently it has been established that the Kondo effect 
plays an important role in transport properties of quan- 
tum dot systems at low temperatures.-^"^ A new feature 
of the Kondo effect in quantum dot systems in compar- 
ison with the traditional magnetic impurity problems is 
that nonequilibrium steady state is realized under a finite 
bias voltage. 

In order to study the properties of the steady states 
theoretically, analyses based on Keldysh formalism are 
often employed. There are two different types of ap- 
proaches to quantum dot out of equilibrium: approaches 
from the noninteracting limit such as the perturbation 
theory with respect to the Coulomb interaction U,^^^ 
and approaches from the strong coupling limit such as 
noncrossing approximation,^ real-time diagrammatic for- 
mulation^*^ and renormalization group method applied to 
the s-d model.^^' 

In the absence of a magnetic field, an anomalous peak 
structure other than the zero bias peak was predicted 
in differential conductance by the 4th order perturba- 
tion theory with respect to the Coulomb interaction UJ 
Recently this new peak was obtained also in the func- 
tional renormalization group result. '^'^ The peak appears 
when the bias voltage exceeds the Kondo temperature, 
eV ~ ksTK, and it is a challenging problem to describe 
the crossover theoretically. 

In finite magnetic fields, it is established both experi- 
mentally and theoretically that the zero bias peak splits 
into two and the peaks are located where the bias voltage 
is equal to the Zeeman splitting, eV = ig/is/i.^-'*'*'^'^'* 

However, the parameter range where each type of the- 
ories can be applied is limited. For example, the per- 
turbation theory is applicable in a relatively weak cou- 
pling regime^ and the noncrossing approximation^ and 
the renormalization group method^^'^^ are confined to 
the strong coupling regime. In the present situation a 
better treatment is required, that covers all the parame- 



ter space of Coulomb interaction, bias voltage and mag- 
netic field. 

Apart from Keldysh formalism, it seems that there are 
only few numerical studies on quantum dot under a fi- 
nite bias voltage. This is largely due to a lack of reliable 
numerical techniques for the nonequilibrium problems of 
the mesoscopic systems. 

One possible numerical approach to the quantum dot 
out of equilibrium is the adaptive time-dependent density 
matrix renormalization group (TdDMRG) method, ^^'^^ 
by which one can accurately calculate time evolution of 
wave function of one dimensional system. In the previous 
TdDMRG study on a quantum dot system^^ the results 
were practically limited in the linear response regime. 
So the V dependence of the physical quantities in the 
nonequilibrium steady state, such as the differential con- 
ductance, were not discussed in their report. 

In this paper we investigate the zero temperature 
transport properties of the quantum dot in a wide range 
of bias votage by the TdDMRG method. We demonstrate 
that, in a limited time interval, the current obtained by 
the finite system calculation provides us accurate infor- 
mation about the one in the steady state of the infinitely 
long system. Then we show that the V dependence of 
the physical quantities can be reliably obtained from the 
time intervals. 

2. Model and Calculation 

8. 1 Hamiltonian 

We consider one single level quantum dot with two 
leads. This system is described by the Anderson model. 
The universal feature of the Kondo effect allows us to 
use 1-D nearest-neighbor tight-binding model for the lead 
parts. In this paper we concentrate on the symmetric case 
for simplicity. Then the model to be studied is described 



H = -t 



i<-l,0<i 



{c]^Ci+ia + h.C.) 
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+ ^{-^ cl^coa + Ucl^co-fcl^coi, (1) 



where the quantum dot is placed at the 0th site. In the 
above equation t is the hopping amphtude in the leads, 
t' the hopping amplitude between the leads and the dot, 
U the Coulomb energy and ha the Zeeman energy. The 
one particle energy at the dot site is set to —U/2 under 
the symmetric condition. 

2.2 Process to reach steady state 

In Kcldysh formalism, one starts with an equilibrium 
state of the unperturbed Hamiltonian. One then turns on 
the perturbation term adiabatically, and gets a nonequi- 
librium steady state after sufficiently long time. 

There are several choices of the perturbation terms 
to be turned on. It is expected that the system reachs 
the same steady state independent of the choices. For 
the present problem, Anderson model out of equilib- 
rium, there are two frequently used options: one is to 
change the hopping amplitudes t' adiabatically, keeping 
the chemical potentials of the two leads constant in time. 
The other is to set t' constant and change the chemical 
potentials. In Kcldysh-bascd theories the former one is 
often used. Note that the interaction term is also turned 
on simultaneously. Here we take the latter choice by 
adding the time-dependent bias term. 



eV 



,i<0 



i>0 



, (2) 



to the Hamiltonian eq.(l). r represents the real time vari- 
able and 0(t) is a smoothed step function which we define 
as 



0{t) = 



exp 



(^) 



+ 1 



(3) 



The reason for our choice of the time dependence of the 
Hamiltonian is that when t' = the number of nonzero 
eigenvalues of the reduced density matrix of the lead 
parts is only one, so the usual DMRG procedure to ob- 
tain the optimal basis set cannot be used directly. 

2.3 Calculation of current 

Based on Keldysh formalism for the quantum dot out 
of equilibrium, the steady current at T = is calculated 
by« 



JiV) 



E 



diV-^^^^^pAco), (4) 

i Lcr + i Her 



where ij-l, l^n arc the chemical potentials of the two leads, 
^L,RiT the resonance widths due to the mixing between 
the dot and the leads, and P(t(w) the spectral function at 
the dot site. For the present model 



(5) 



where Dl^r„{lo) is the local density of states at the edge 
of the lead. The bias voltage V appears in T l,r, Pa{'jj) 



and as hl — PR = eV in eq.(4). The essential part to 
compute eq.(4) is reduced to the calculation of P(r(a;). 

On the other hand, in the TdDMRG calculation we 
simply get the current by taking the expectation value 
of the current operator with the wave function at each 
time. We may define two current operators: one flows 
from the left lead to the dot and the other from the dot 
to the right lead, 

Mr) = (V>(r)|i^5^(cLc-i. -/i.c.)|V(t)), (6) 



2h 

a 

Mr) = {Hr)\'-^Y.(^lco. - h.c.mr)). 



(7) 



In the TdDMRG the essential part to obtain the current 
becomes the calculation of \iI){t)). This can be done by 
the TdDMRG method. 

In the noninteracting case, we can easily compute 
Po-(a;) exactly and obtain J{V) for infinitely long leads. 
Alternatively, when U = we can also use exact di- 
agonalization for the Hamiltonians without the bias 
term and calculate \iP{t)) exactly by assuming sudden 
switching-on of the bias voltage. We use these results to 
compare with the TdDMRG results in later sections. 

3. Adaptive Time-dependent DMRG Method 

Here we briefly explain the numerical technique we use. 
The basic idea is to combine the Suziiki- Trotter decom- 
position of the time evolution operator and the DMRG 
finite-system algorithm with the wave function predic- 
tion method. "'^^ 

Our Hamiltonian eqs.(l) and (2) can be written as the 
sum of the ith bond Hamiltonian hi{T). Then the time 
evolution operator is found to be 



U{to + At,to) 



To+Ar 



hi{T)dT 



~ e 



X e ^ 2 



(8) 



by using the 2nd order Suzuki- Trotter decomposition of 
the T-ordered exponential.^^ Then the time evolution op- 
erator is well approximated as the product of the local 
time evolution operators [7,!°'^''*'. 

The ground state wave function obtained from the 
usual DMRG calculation for open boundary condition 
has the form of 



E 



(9) 

where \ai), \ai+\), |cr/+2), lA-i-s) are the DMRG basis set 
of the left block (which consists of I sites), left site, right 
site and right block, respectively. The operation of the 
local time evolution operator at the sites I + 1 and I -\- 2 
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can be calculated without any error as 

V l + i- ^/Ql<Ti + lCTi_|_2Pi + 3 

(10) 



E 



/ T /-local \ 



Then the wave function prediction method is used to 
move to the next configuration of the finite system algo- 
rithm, that is, the number of the basis to describe the 
system block and the one additional site is truncated to 
a fixed number m. 

Performing this procedure at every step of the DMRG 
finite system algorithm, the full operation of the decom- 
posed time evolution operator eq.(8) on the wave func- 
tion can be done in one full sweep, keeping the basis set 
optimal. 

The main sources of errors involved in this method are 
the Trotter error, which comes from the Suzuki- Trotter 
decomposition of the time evolution operator, and the 
DMRG truncation error, which originates from the re- 
duction of the number of block basis to m in the wave 
function prediction method. As the calculation proceeds 
the errors acummulate step by step and the TdDMRG 
result gradually loses its accuracy. 

The TdDMRG results shown in this paper are ob- 
tained with At — 0.05h/t. 

4. Time Dependence of the Current 

Since we are interested in the properties of the 
nonequilibrium steady states of the infinitely long sys- 
tems, the most important requirement for the present 
TdDMRG calculation is to realize the steady states nu- 
merically. But it is obvious that the steady states cannot 
be realized in a finite system in a rigorous sense. 

In this section we show that one can find steady- 
like behaviors in finite systems that mimic behaviors of 
steady states of infinite systems and discuss how to ob- 
tain the information on the steady states from the Td- 
DMRG calculations. 

4-1 Oscillations in Jl{t) and Jr{t) 

The total particle number and are the conserved 

quantities of the Hamiltonian eq(l) and eq.(2). Thus it is 
sufficient to consider only the states which belong to the 
subspace with fixed particle number and J2i ■ ^tii^ 
paper we concentrate on the half-filled and 5f — 
subspace, where the Hamiltonian is the most symmetric 
and the Kondo effect is the most significant. Then the 
system length L should be taken as an even number and 
the lengths of the left and right leads cannot be the same. 
We set the right lead is longer than the left lead by one 
site. 

Jl{t) and Jr{t) defined in eqs.(6) and (7) should be 
the same in steady states in infinite system. However for 
a finite system these currents show oscillations^" as in 
Fig.l. It can be seen that the amplitude of the oscillation 
is smaller for L = 400 than for L ^ 64 and therefore 
we conclude that this oscillation disappears in the limit 
L — > oo. For the present model under the symmetric 




Fig. 1. (Color online) J/, (t) and Jfl(T) with different system sizes 
(L = 64 and L = 400) obtained by the exact diagonalization. The 
parameters used are t' /t = 0.6, U/t = 0, eV/t = 1.0 and h/t = 0. 
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Fig. 2. (Color online) J(t) for L = 64 and L = 400 obtained 
by the exact diagonalization. The parameters used are L = 
400, t'/t = 0.6, U/t = 0, eV/t = 1.0 and h/t = 0. 



condition Jl{t) and Jr{t) oscillate in opposite phase. 
Consequently we use the averaged current 

1 



Jir) = 7;{Jl{t) + Jr{t))^ 



in order to remove the oscillations. 



(11) 



4-2 Behavior of the current in long time scale 

The long time behavior of J(t) is dominated mainly by 
the system size L. In Fig. 2 we see an obvious difference 
between L = 64 and L = 400 results after r ~ 30h/t. 
This can be explained as follows: electron wave packets 
driven by the bias voltage will eventually arrive at the 
right edge of the system, and be reflected without any 
loss of energy. After some time the wave packets will get 
back to the center of the system and make negative con- 
tributions to J(t). J(t) for L = 64 after r ~ 30h/t re- 
sults from a complicated superposition of the wave pack- 
ets going left and right. 

Because this reflection at the edges of the system is an 
artifact of the flnite size calculation, we do not consider 
J(r) after this effect appears. 

4-3 Switching-on of bias voltage and response 

A characteristic energy of the system in equilibrium is 
the rcnormalized resonance width F = F/x|-f, where 
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is defined as^^ 



Xo 



d{ha') 



(12) 



using the retarded self energy for V — Q, Xaa' in 

equilibrium state can be calculated by the Bethe Ansatz 
solution, r is a monotonically decreasing function of U 
and becomes the Kondo temperature Tk times a con- 
stant in the strong coupling limit. 

Response of the system to the switching-on of bias 
term takes place in the time scale of the order of ?i/r. 
When the time variation of the bias term is faster com- 
pared to fi./r, J(r) shows an overshoot behavior fol- 
lowed by a damping oscillation. Therefore in order to get 
smoother results for J(t), it is preferable to use slower 
switching-on of the bias term. On the other hand, the 
result in long time scale becomes unreliable because of 
the reflections at the edge of the system and the acum- 
mulation of errors. 

In an actual experimental situation the resonance 
width r is much smaller than the band width 4i. However 
we take rather large value of F, with the aim to acceler- 
ate the response of the system and to lead the system to 
reach the steady state before the calculation loses its re- 
liability. We fix t' jt as 0.6, which means r(w — 0)|y=o — 
0.72t. In addition we set tq = 3H/t,Ti — h/t. We will 
see later that we can obtain reasonable results by using 
these values of t' and the parameters in the smoothed 
step function. 

4-.^ Quasi- steady state in finite system 

As already stated, in the noninteracting case we can 
use three different methods to calculate current: inte- 
gration of eq.(4), exact diagonalization and TdDMRG 
method. In Fig. 3 we compare the results obtained from 
these methods. For U /t = the three results show 
good agreement in the flat region: Sh/t < t < 30h/t 
for eV/t = 0.5, 8h/t < t < 2?,h/t for eV/t = 1.0, 
Sh/t <T < im/t for eV/t = 1.5 and 8h/t < t < Ibh/t 
for eV/t = 2.0. Note that the system sizes and the 
processes to bring the system to the steady state are 
quite different. Nevertheless the currents in the flat re- 
gion show a good coincidence. Thus we can conclude that 
the system reaches similar steady state independent of 
the choice of the perturbation term. Moreover, because 
L = oo for the analytic calculation and L — 64 for the 
TdDMRG, we may also conclude that the system size de- 
pendence of the steady current is small under the sym- 
metric condition. Out of the symmetric condition this 
system size dependence is considerably large, and one 
has to employ the finite size scaling to obtain results for 
the infinite system. 

Even for U/t = 1 or 2 we can get the well defined 
plateaus of J(t) in Fig.3. For U/t = 2, eV/t ^ 2 the 
TdDMRG result with m = 600 does not stay at a par- 
ticular value but with m = 800 the plateau appears from 
r ~ 9h/t to T ~ lAh/t. In these plateau regions the sys- 
tem is expected to simulate the nonequilibrium steady 
states of the infinitely long system, just as in the U = 
case. 
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Fig. 3. (Color online) J(t) obtained by analytic calculation (L = 
oo.U = 0), exact diagonalization (L = 400, f/ = 0) and Td- 
DMRG (L = 64, C//t = 0,1,2). The steady-like behaviors can 
be seen in all the TdDMRG results shown here except for the 
one U/t = 2, eV/t = 2, m = 600. The parameters used are 
t'/t = OS, h/t = and eV/t = 0.5,1.0,1.5,2.0 from the top. 



Hereafter we call the system in the fiat region where 
J{t) stays almost constant within ±1% as a quasi-steady 
state in finite system. Then it is natural to define the 
steady current as an average value over the quasi-steady 
state region. We will see later that this definition gives 
consistent results with the Keldysh results for not only 
U — but also J7 7^ at least in the low bias regime, see 
§5.1.1. 

It is worth noting that the overshoots discussed in the 
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previous subsection are observed in Fig. 3. We see that 
the overshoots become bigger for larger U and V. This 
is due to the small value of F, the energy scale of the 
system, relative to the time variation of the bias term. 
These overshoots are expected to disappear if we use 
more slowly-changing bias voltage, in other words, if we 
use larger values of tq, ti. But in this paper we do not try 
to fine-tune the optimal choice to reduce the overshoots, 
because the quasi-steady states are well-defined in the 
parameter range we have studied. 

4-.5 Decay of quasi-steady state 

As we have seen, we cannot realize the true steady 
states by the TdDMRG method. We can only realize the 
quasi-steady states that have finite lifetimes. 

There are several factors that determine the end of the 
quasi-steady state. One is the refiection of the electron 
wave packets at the edges of the system, as discussed in 
§4.2. Another is due to the change of the population of 
electrons in each lead. When the number of electrons in 
the left lead becomes too few, the current cannot keep 
on flowing steadily. The third one is the acummulation 
of the truncation error, see §3. 

The larger bias voltage V means the more current flows 
and the larger change occurs in the wave function by 
the time evolution. Then the second and third factors 
become important when V is large. Therefore, the larger 
V we use, the shorter the lifetime of the quasi-steady 
state becomes. From Fig. 3 we can confirm that the decay 
happens earlier with increasing bias voltage. Moreover we 
can see that by increasing m, in other words by reducing 
the truncation error, the decay of the quasi-steady state 
can be delayed as in the results for U/t = l,eV/t = 1 
and U/t = 2, eV/t = 2. 

Note also that from the results U/t — 1, eV/t = 1, m = 
600 and m = 800, the m dependence of the current in 
the plateau region is small. Thus we can expect that once 
the clearly recognizable quasi-steady state is realized the 
steady current obtained is very close to the exact value 
for m oo. 

Eventually it is seen that the quasi-steady states are 
well realized for <U/t < 2, 0.5 < eV/t < 2 by using 
the parameters chosen. 

4-. 6 Boundary conditions 

It is well known that the Kondo cloud spreads as U 
increases. In other words, the lower energy excitations of 
the leads become important with increasing U. 

From this reason one has to employ sufficiently large 
system to describe the Kondo effect successfully by 
DMRG. However it is numerically difficult because a 
large system requires a large number of basis to be kept 
for a given accuracy. Then the calculation for larger U 
takes much more time. 

For the present calculation, the low bias regime is dif- 
ficult to treat because the low energy excitations are 
mainly involved in the transport, and thus the size effect 
is strong. We can improve effeciency of the calculation 
in the parameter range by using a different boundary 
condition other than the usual open boundary condition 
(OBC).^^' The idea is to make the energy scale near 
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Fig. 4. (Color online) TdDMRG results for J(t) with different 
boundary conditions. It can be seen from the OBC result that 
the system size is not sufficient to describe the Kondo effect. 
Obviously the SBC result is much better than the OBC one in 
the sense of finding the quasi-steady state. The parameters used 
are L = 64, m = 600, t'/t = OS, U/t = 2.0,eV/t = 0.01 and 
h/t = 0. Inset: The hopping amplitudes with SBC near the left 
edge of the system obtained by eq. (13). 



the boundary of the lead small by reducing the hopping 
amplitude exponentially towards the edges of the system, 
thus including the lower energy excitations. There are 
of course several possible ways of reducing the hopping 
parameters. In this paper we use the smooth boundary 
condition^^ (SBC) by setting the 10 hopping parameters 
from the edge as 

t.,,,A^-(^l + tanh-^j (13) 

where Xi ^ {i — l/2)/10, as shown in the inset of Fig. 4 
(left side). In Fig. 4 we see dramatic improvement in the 
realization of the quasi-steady state. While the OBC re- 
sult does not show the steady-like behavior, we clearly 
observe a quasi-steady state for the SBC result. This im- 
phes that our choice of the hopping parameters (eq.(13)) 
works well for the parameters in Fig. 4. 

It should be mentioned that there are some short- 
comings of using the SBC. First is that the SBC slows 
down the convergence of the diagonalization step in the 
usual DMRG method. This is due to the facts that the 
diagonalization of a large sparse matrix is done by a 
power method such as Lanczos or Davidson method, and 
that the SBC bring the first excited energy closer to the 
ground state energy. The second is that in the SBC re- 
sults the truncation errors become more significant than 
in the OBC resuhs. So we use the SBC for eV/t < 0.3, 
otherwise we use the OBC. By setting the boundary con- 
ditions in this way, we can obtain the quasi-steady states 
for 0<U/t<2,0< eV/t < 2. 

4-7 1/L correction for h ^ 

Thus far we have concentrated on the zero magnetic 
field case, where the system size dependence of the phys- 
ical quantities is small due to the symmetric condition. 
By making h finite, small corrections appear in the val- 
ues of steady current J{V) and expectation values of 
n-f = cJ|CoT both in equilibrium n-|-(T = 0) and in the 
quasi-steady state ni{V) as is seen in Fig. 5 for the non- 
interacting case. Here the oscillation in ?t.|(t) comes from 
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Fig. 5. (Color online) J(t) (upper figure) and n|(T) (lower figure) 
with different sistem sizes obtained by the exact diagonalization 
(L = 64, L = 400) and analytic calculation (L = oo). The pa- 
rameters used are t' /t = 0.6, U/t = 0, eV/t = 1.0 and h/t = 0.5. 



Fig. 6. (Color online) The current obtained via eq.(4) for L = 
OD,U/t = and the current calculated as the averages over the 
quasi-steady state region of the TdDMRG results. The param- 
eters used are L = 64, h/t = Q, m = 800 for U/t = 2, eV/t = 2 
and m = 600 for the others. 




the oscillations in Jl{t) and Jr{t), and we define n]{V) 
as the central value of the oscillation in the quasi-steady 
state region. The finite size correction in J{V) has a plus 
sign and the correction in n^{V) has a minus sign. The 
system size dependence of these corrections in J{V) and 
n-^iV) is found to be 1/L. 

This phenomena can be explained as follows: when the 
spin polarization at the dot (S^) becomes finite by the 
symmetry breaking caused by the magnetic field, the po- 
larization of the lead part as a whole is — (S*^), since 
J2i Si is fixed to zero. Because of this finite spin polar- 
ization the lead part of a finite system is shifted from the 
siglet state, which is the ground state of the lead part, 
whereas for an infinite system —{S^) does not cause any 
effect to the lead part. {S^) of a finite system is deter- 
mined energetically by the balance between the Zeeman 
energy and the energy difference of the lead part from the 
siglet state. The latter is important for relatively short 
system, suppressing (5^). Thus, it can be concluded that 
the magnetic field is effectively reduced by the finite size 
effect. This effect makes J{V) larger compared to that 
of longer system. 

In the interacting case, it is expected that this finite 
size effect becomes stronger than the U = case because 
of the enhancement of the polarization by the Coulomb 
repulsion. 

Though we can remove this l/L correction by the finite 
size scaling, we simply neglect it in this paper. Hence 
the corrections are included in the h ^ results in later 
discussions, but this is not important for < U/t < 2. 



Fig. 7. (Color online) The TdDMRG current and J(V) analyt- 
ically calculated with the use of eqs.(4) and (14), both for 
U/t = 2. Note that the latter is asymptotically correct in the 
low bias regime. The TdDMRG result here is the same data set 
shown in Fig. 6. 



5. Physical Properties of the Steady State 

In this section we discuss the dependence of physical 
quantities on the Coulomb interaction U, the bias voltage 
V and the magnetic field h. The results presented in this 
section are obtained as follows: for each value of V we 
perform the TdDMRG calculation, get the quasi-steady 
state and take the average value over the quasi-steady 
state region. 

5. 1 Zero magnetic field 
5.1.1 Current 

Fig. 6 summarizes our results of JiV). 

In the noninteracting case we see nice agreement be- 
tween the TdDMRG and the analytic results obtained 
via eq.(4), for all bias voltages. 

For t/ / let us compare the TdDMRG results with 
analytic ones in the low bias regime. For this purpose we 
use an asymptotic expression for p in the low energy, low 
bias voltage and low temperature,^** 




Putting this expression into eq.(4), we obtain the asymp- 
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Fig. 8. (Color online) The differential conductance for various 
U/t obtained by the integration of eq. (4) (U/t = 0) and the 
TdDMRG and interpolation. The parameters used are the same 
as in Fig. 6. 



totic behavior of J(V) for small V. From the compari- 
son for U/t = 2, Fig.7, we see that the TdDMRG re- 
sult is in excellent accordance with the analytic one at 
eV/t < 0.25. Therefore we can conclude that our Td- 
DMRG calculation and definition of the quasi-steady 
states work well for U/t < 2 at least for low V. 

Back to Fig. 6, in the high bias voltage we see J{V = 
1.75t/e) > J{V = 2t/e). This behavior is due to the 
energy dependence of T^ jf , thus representing a specific 
nature of the 1-D nearest-neighbor tight-binding leads. 
Since we are interested in the universal behavior of the 
Kondo effect in the quantum dot system, we will restrict 
ourselves to study up to eV/t = 2. 

5.1.2 Differential conductance 

In many works on the nonequilibrium transport phe- 
nomena, the differential conductance G{V) = ^^y-^ has 
been discussed as an extention of the linear conductance. 

Here we carry out the following process to get G{V) 
from the data points of J{V). We interpolate the J{V) 
for < eV/t < 1.5 by the least square fitting to an odd 
polynomial expression. Then we differentiate the result- 
ing function and obtain G{V), Fig. 8. Note that G{V) 
calculated in this way contains errors due to the inter- 
polation especially near eV/t = 1.5, the edge part of the 
data points. 

In Fig. 8 the unitarity limit of the zero bias peak 
G{0) — 2e'^/h are clearly seen for every U. G{V) drops 
from the zero bias peak with increasing V. The width 
of the zero bias peak becomes narrower as U increases. 
These behaviors result from the well known fact that the 
width of the Kondo peak in Pa{^) is essentially given by 
the Kondo temperature Tk, and that the peak height of 
the Kondo peak is suppressed by the bias voltage. Fig. 8 
demonstrates that the TdDMRG succesfully reproduces 
the characteristic properties of G{V).^''^'^ 

5. 2 Effect of magnetic field 

5.2.1 Current and differential conductance 

Even under a finite magnetic field up to < 0.5 our 
TdDMRG calculation does not lose its validity to study 
the steady state, in spite of the 1/L correction discussed 
in §4.7. Using the same definition of the quasi-steady 
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Fig. 9. (Color online) The current under a magnetic field for 
U/t = 0,1,2 from the top. The top {U/t = 0) result is ob- 
tained by the integration of eq. (4) and the others are obtained 
by the TdDMRG calculation. The results for h/t = are the 
same ones in Fig. 6. The other TdDMRG results are obtained 
with L = 64, m = 1024 and the OBC. 



state and the same interpolation as in the h = case, we 
get the resuhs for J{V), Fig.9 and G{V), Fig.lO. 

It is seen that J{V) is suppressed by the magnetic 
field for all U and V. Then we observe the splitting of 
the zero bias peak in G{V) by a finite magnetic field 
when U and h are relatively large. Moreover the position 
of the splitted peak is roughly equal to the twice of the 
Zeeman energy, eV = 2h. These behaviors are consistent 
with the previous results,^' ^ and can be naively explained 
by the width of the peak of Pai^) and the separation 
of p^{uj) and pi{uj) by the Zeeman energy. When U is 
large the peak is narrow, then the separated peak is not 
included in the interval of integration in eq.(4), if eV 
is small compared to the Zeeman energy. Thus in this 
case G{V) has the splitted peaks corresponding to the 
peaks of Pcr{'^)- On the other hand, when U is small the 
separation of the broad peaks in P(j{^) does not affect 
the peak structure of G{V). 

Again we can conclude that our TdDMRG calculation 
works properly for the calculation of G{V) in magnetic 
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Fig. 10. (Color online) The differential conductance under a mag- 
netic field for U/t = 0, 1, 2 from the top, obtained from the data 
shown in Fig. 9. 



Fig. 11. (Color online) The magnetic susceptibility defined as in 
eq.(15) for U/t = 0,1,2 from the top. The figure at the top 
[U/t = 0) is obtained by eq. (16). The others are obtained by the 
TdDMRG with L = 64, m = 1024 and the OBC. 



fields. 

5.2.2 Magnetic properties 

With the knowledge of [■(/'(t)}, of course, one can take 
expectation values of various operators other than the 
current operator. In this paper we present here the sus- 
ceptibility, Fig.ll, which we define as 

xiV) ^ (15) 
h 

where magnetization at the dot (S^) is the averaged value 
of (V'(t)|5^|?/'(t)) over the quasi-steady state region. The 
numerically calculated {ip{T)\ni\i/j{T)) shows the oscilla- 
tion as in the lower figure of Fig. 5, but the oscillation can 
be removed by subtracting {^p{T)\ni\'ip{T)). In Keldysh 
formalism the magnetization can be calculated from 

/OO J 
(16) 

where G~~^{uj) is the Fourier component of one of the 
Keldysh Green functions, G^+(t) = i(cJ^(0)coo-(T)). 

Let us start the interpretation of Fig.ll from V = 0. 
It is well known from the analyses of the Kondo effect 



in equilibrium that the Coulomb repulsion U enhances 
the linear susceptibility lim^^^ x(^ = 0) and a finite 
magnetic field suppresses x{V = 0). We can clearly see 
these behaviors in Fig.ll. 

Next the effect of V is explained qualitatively by con- 
sidering the simplest case, where U — and the density 
of states of the leads is constant. Then we easily find 

{S^)^l J '''''\u{p^{uj)-p^{u)), (17) 

using the symmetric condition. From this expression we 
see that when h ^ eV the peak of p| is included in the 
integration but not the peak of px, resulting in a large 
value of (S*^). On the other hand, p| and pi are both 
small in the interval of integration when h <^ eV . Thus 
we can expect that large V suppresses the spin polar- 
ization and the susceptibility, as can be seen in Fig.ll. 
Even when U ^ and the density of states has energy 
dependence, the above story roughly holds and explains 
the V dependence of x(F) in Fig.ll. 



J. Phys. Soc. Jpn. 



Full Paper 



9 



5.3 Behavior of (n^ni) 

In this subsection let us discuss the U, V and h depen- 
dence of (n-^n^). Since this operator appears as Un-^ni 
in the Hamiltonian (1), its expectation value reflects the 
effect of electron correlation. 

Again (?/'(T)|nf nJ'0(T)) as a function of time shows 
an oscillation as in Jl.r.{t) and n^^ilr). In this case 
we remove the oscillation by taking the average value 
of (Sn-^Sni) over the quasi-steady state region, where 
6na —ricr - {ria). 

From the previous study it is known that in the high 
voltage limit the third-order contributions to the self en- 
ergies vanish, resulting in 1/2.'* This is accounted 
for by effective suppression of electron correlation and 
magnetic field by the bias voltage. Thus it is likely that 
— * 1/4 in the limit of I^ ^ cxd, independent of h. 
In Fig. 12 we see for U/t = 0, approaches to 1/4 
monotonically with increasing V . However, for U/t =1,2 
there appears nonmonotonic behavior for relatively small 
/i, while the magnetization is a monotonically decreasing 
function of V , as can be seen in Fig.ll. In other words, 
{n-^ni) has a minimal value at eV/t ^ 1.2. Consequently 
the electron correlation seems to be enhanced effectively 
by the bias voltage, compared to the equilibrium state. 
This behavior is new and needs an explanation. 

This anomalous behavior is seen in relatively low volt- 
age regime eV/ 1 < 1 and the effect of the energy depen- 
dence of DL,B.a, the local density of states of the leads, 
becomes important in the higher bias voltage regime. 
Therefore we may expect that this new behavior is a uni- 
versal behavior of the Kondo effect of the noncquilibrium 
steady state in the quantum dot system. 

6. Conclusions 

In this paper the adaptive time-dependent DMRG 
method was applied to the 1-D Anderson model with 
time dependent bias term. We have shown that the Td- 
DMRG works well for the studies of the quantum dot 
out of equilibrium. 

In §4 it was seen that from J(t) for < U/t < 2, < 
eV/t <2, < h/t < 0.5 obtained by the TdDMRG 
calculation the quasi-steady states were found. In order 
to realize the well-defined quasi-steady states, we made 
some technical improvements compared to the previous 
study:^^ first, we simply made the number of the DMRG 
basis m larger to keep the accuracy until sufficiently long 
quasi-steady state appeared in J(t). Second, we used 
the slowly changing bias voltage other than the suddenly 
changing one, and reduced the effects of the overshoots 
and the damping oscillations in J(r). Third, for low bias 
regime the SBC was used to obtain flat J(t). By these 
improvements, the studies of the properties of the quasi- 
steady states became possible. 

Then we verified that the physical quantities averaged 
over the quasi-steady state regions behave exactly like 
the ones of the real noncquilibrium steady states in the 
infinite system, except for the 1/L corrections for h ^ 0. 
The TdDMRG calculation successfully reproduced the 
established results for the differential conductance, such 
as the zero bias peak G(0) = 2e^/h for h = 0, and the 
splitting of the zero bias peak caused by a magnetic field 
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Fig. 12. (Color online) The expectation values of n^n^ taken by 
the quasi-steady state wave function for U/t = 0, 1, 2 from the 
top. The top figure [U/t = 0) is obtained by (16) and {n-^n^) = 
(n|)(n|). The others are obtained by the TdDMRG calculation 
and the parameters used are the same ones as Fig.ll. 



with the peak position eV — 2h. Also for x(^) ob- 
tained reasonable results. 

Additionally we observed a nonmonotonic behavior in 
(n-fTij^) as a function of V. We discussed that this may 
be a universal behavior of the Kondo effect in quantum 
dot system. 

It is concluded that the TdDMRG method is a useful 
numerical tool to investigate the noncquilibrium steady 
states of quantum dot systems. However it should be 
noted that there is a restriction to this method due to 
the large Kondo cloud in the strong coupling regime. 
In order to correctly describe the Kondo physics in the 
strong coupling regime, one must employ a very long 
system and therefore one has to take a very large value 
of m. In this paper we studied relatively weak coupling 
regime U/t < 2 and use the SBC to control this size ef- 
fect. For larger U it is expected that much more time is 
required to correctly realize the quasi-steady states espe- 
cially in the low bias regime. Except for this limitation, 
the method can be applied to other nonequilibrium prob- 
lems of mesoscopic systems. 
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